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ABSTRACT 

The model of magne tic braking of solar rotation considered by 



Charbonneau fc MacGregorl (119931 ) has been modified so that it is able to re- 
produce for the first time the rotational evolution of both the fastest and slow- 
est rotators among solar-type stars in open clusters of different ages, without 
coming into conflict with other observational constraints, such as the time evo- 
lution of the atmospheric Li abundance in solar twins and the thinness of the 
solar tachocline. This new model assumes that rotation-driven turbulent diffu- 
sion, which is thought to amplify the viscosity and magnetic diffusivity in stellar 
radiative zones, is strongly anisotropic with the horizontal components of the 
transport coefficients strongly dominating over those in the vertical direction. 
Also taken into account is the poloidal field decay that helps to confine the width 
of the tachocline at the solar age. The model's properties are investigated by 
numerically solving the azimuthal components of the coupled momentum and 
magnetic induction equations in two dimensions using a finite element method. 

Subject headings: MHD — stars: evolution — stars: interiors — stars: magnetic 
fields — Sun: rotation 



Introduction 



Helioseismology has revealed that the Sun's radiative core rotates nearly uniformly, 
at least down to the radius r ~ 0.2 Rr? u with a rate that is close to th e mean rotation 



rate of its convective envelope (e.g., iTomczyk. Schou. fc Thompson! Il995t ICouvidat et al. 



20031 ). It is also known that the Sun, like other solar-type stars, has bee n experiencing an 



effective braking of its surface rotation via a magnetized solar (stellar) wind ( )Schatzmanlll962 
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Weber fc Davisl 119671 ; iKawalerl Il988l ). Because the magnetic field lines that sling charged 
particles into space are rooted in the photosphere, it is the envelope rotation that should be 
decelerated by the wind torque, while the conservation of angular momentum in the radiative 
core should keep it rapidly rotating. This expectation appears to be in a stark contrast with 
the measurements of internal solar rotation, which indicate that the difference in the angular 
velocity between the Sun's radiative core and convective envelope must have been erased by 
some angular momentum redistribution mechanism. A similar conclusion can be reached 
from comparisons of the rotation period distributions for solar-type stars in open clusters of 
different ages. It turns out that the core and envelope rotation have to be coupled with a 
characteristic timescale of angular momentum transfer between them inc reasing from a few 



Myr to n early a hundred Myr between the fastest and slowest rotators (IDenissenkov et al. 



2009 


; cf. 


Irwin et al. 


2007) 



angular momentum transport described as a purely diffusive process shows that this coupling 
time interval roughly corresponds to a diffusion coefficient decreasing from 10 6 cm 2 s _1 to 
5 x 10 4 cm 2 s _1 . However, such large diffusion coefficients would lead to rates of surface Li 
destruction strongly exceeding the trend observed in the Sun and its twins (Section 12. 4p . 
Therefore, we have to admit that the angular momentum redistribution in solar-type stars is 
unlikely to be accompanied by an equally rapid mass transfer. This restricts our search for 
possible mechanisms of angular momentum transport in solar-type stars to those employing 
waves and magnetic fields. 

Presently, the most popular mechanisms of angular momentum redistribution in the 
Sun's radiative core are the smoothing of its rotation profile b y internal gravity waves (g- 
modes) generated by turbulent eddies in the c onvective envelope (ICharbonnel fc Talonll2 005'l 



and magne tic braking of differential rotation ( jCharbonneau fc MacGregor 



1992 



1993 



Spruit 



19991 . 120021 ). Using a simplified a nalytical prescription f or the spectrum of internal gravity 
waves and other approximations, ICharbonnel fc Talonl (120051 ) have succeeded in obtaining 
a nearly flat rotation profile in the present-day Sun simultaneously with a time evolution 
of its internal rotation consistent with the observed solar Li abundance. However, because 
of the complexity of physical processes involved in the generation, propagation and dis- 
sipation of internal gravity waves in the Sun, it is difficult to assess the validity of this 
solution. For instanc e , it w ould not be valid if a flatter spectrum, like that computed by 
Rogers fc Glatzmaierl ( 120051 ). were used. Non-linear wave- wave interactions, as well as the 
i nteraction of gravity waves with a magnetic field, could also change t he solution drastically 



( jRogers. MacGregor. fc Glatzmaier 



2008 



Rog ers fc MacGregorl 120091 ) . Even in a case very 



similar to that considered by ICharbonnel fc Talonl ( 120051 ) the question arises as to why grav- 
ity waves do not disturb the uniform solar rotation? Indeed, the power of gravity waves 
generated in the present-day Sun should approximately be the same as it was in the young 
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Sun; consequ ently, given their peculiar anti-diffusive nature in shear flows (e.g.. lPlumblll977 
Ringotlll9981). gravity waves should have forced the solar r otation profile away from the uni- 
form one (IDenissenkov. Pinsonneault. &: MacGregorl 120081 ) . However, this has not occurred. 



For a magnetic field configuration to influence differential rotation of the solar radiative 
core, it must have a non- vanishing azimuthal (toroidal) component. To achieve this, it is 
usually assumed that there is an axially symmetric poloidal magnetic field frozen into the 
plasma inside the radiative core which can be stretched around the rotation axis by the 
differential rotation itself to form the necessary toroidal field. As far as the origin of the 
poloidal field is concerned, it is widely believed that the field was either inherited from a 
protostellar cloud or generated by a convective dynamo during the Sun's pre -MS evolution 
when its convective envelope occupied a much larger volume. ISpruitl ( 119991 ) proposed the 
original alternative explanation that the poloidal field might be constantly replenished by 
radially stretching the toroidal field in a dynamo-like cyclic process, the necessary radial 
plasma displacements being caused by the non-axisymmetric (m = 1) kink instability of the 
toroidal field (the Tayler-Spruit dynamo). 

Magnetic braking is implemented by the azimuthal component of the magnetic tension 
term of the Lorentz force. Given that its corresponding engagement time is of order several 
Alfven timescales for the po loidal field, magnetic breaking appears to be very efficient even 
for a small field amplitude (IMestel fc Weisslll987l ). The only problem associated with this 



mechanism is that it leaves behind a torus-shaped "dead" zone inside the radiative core 
with the angular velocity increasing toward the center of the torus cross-section. This is a 
direct consequence of Ferraro's law of isorot ation (rotation that keeps the angular velocity 
constant along each field line) ( Ferraro 1937}) and the presence of a constant dipole poloidal 
field. To solve this problem, ICharbonneau fc MacGregorl (11993) (he reafter, referred to as 
CHMG93) and other researchers (e.g., iRuediger fc Kitchatinovl Il996l ) artificially amplified 
the viscous transport with the justification that rotation and associated hydrodynamic in- 
stabilities should lead to turbulence, and hence that the microscopic viscosity should be 
augmented by a much stronger turbulent viscosity. However, the minimum value of the 
enhanced viscosity (5 x 10 4 cm 2 s -1 ) with which it is possible to get a nearly flat rotation 
profile in the present-day Sun results in too much Li being transported below the bottom of 
convective envelope and thus destroyed, inconsistent with the observed solar Li abundance. 
Therefore, the proposed isotropic viscosity enhancement cannot be considered a satisfactory 
solution. 

Only recently has the physics-based modeling of angular momentum transport in solar- 
type stars begun to use the extensive data sets of rotation periods for solar analogs in open 
clusters as an observational constraint. The main goal of previous studies was to demonstrate 
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that a part icular model could (or co uld not) produce the solid-body rotation of the Sun. In 
particular, Eggenberger et al.l (120051 ) showed that the Tayler-Spruit dynamo could account 
for the flat rotation profile of the Sun. However, the correct model must also explain the 
cause of the transition from the solid-body rotational evolution of the fastest rotators to 
the evolution wherein a radial differential rotation is sustained for a hundred Myr in the 
slowest rotators. It turns out that the Tayler-Spruit dynamo cannot pass this observational 
test because it alw ays enforces a nearly uni form rotation in solar-type stars, no matter how 
slowly they rotate ( iDenissenkov et al.ll2009l ). 



Another important problem that the correct mechanism of angular momentum transport 
in the Sun has to explain is the thinness of the solar tachocline. The tachocline is a layer 
in which the differential rotation (both radial and latitudinal) of the con yective envelope 
sharp ly changes to the nearly solid-body rotation of the radiative core f lSpiegel fc Zahn 
1992 ). Its thickness, as measured by helioseismology, is only about 4% of the solar radius 
( jfjasu fc Antia 2003 ). It has been suggested that no purely hydrodynamic mechanism can 
explain its existence and that a large-sca le magnetic field in the radiative core is therefore 
neede d to confine the tachocline structure (jGough fc McIntyrdll998l : lRuediger &: Kitchatinov 
2007h . 



To summarize: because large-scale magnetic fields have certainly been playing crucial 
roles during different phases of the Sun's rotational evolution — such as the synchronizing of 
rotation of the young Sun and its sur rounding protoplanetary disk (the disk locking) (e.g., 
Shu et al.lll994j ; iMatt fc Pudritzll2005l ). providing the leverage for the solar wind, confining 
the tachocline, being generated through a convective dynamo and then driving the solar 
activity — we believe that it is worthwhile to develop a mechanism of magnetic breaking of 
differential rotation of the Sun's radiative core. 

In this paper, we ela borate upon and exten d the magnetic braking mechanism that 
was originally propose d by iMestel &: Weis s (1987) and then studied in detail numerically by 
CHMG93 followed by lRuediger &: Kitchatinovl (119961 ). Our only modification of the mecha- 
nism is to assume that the turbulence induced by rot at ion- driven hydrodynamic instabilities 
in the Sun's radiative core is highly anisotropic with its corresponding horizontal component 
of turbulent viscosity strongly dominating over that in the vertical direction u v . The 
hypothesis o f anisotropic turbulent diffusion in stellar radiative zones was first advanced by 
Zahnl (119921 ) and it had since been productively used by many other researchers to study 



rotat i onal mixing in both MS and red g i ant b ranch stars ( jTalon fc Zahnl 119971 ; iTalon et al. 



19971 : iMaedeilliooi IPalacios et alll2003[ 200 6J). Initially it w as even considered to explain 



the thinness of the solar tachocline by ISpiegel fc Zahnl (119921 ) . We will show that this hy- 
pothesis alone permits a solution of almost all of the aforementioned problems associated 
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with this particular mechanism of magnetic braking, such as the Li problem (a sufficiently 
small value of the diffusion coefficient v y can be chosen so that it does not lead to excessive 
Li destruction), the dead- zone problem, and the difference in spin-down of the fastest and 
slowest rotators among solar-type stars in open clusters. 

The paper is organized as follows. In Section [2j we briefly discuss simple models of 
rotational evolution of solar-type stars and review the main results that were obtained by 
comparison of the model predictions with the observed period distributions for solar counter- 
parts in open clusters. In Section [5TT| we reproduce the 2D MHD computations of CHMG93 
and show that they cannot account for the spin-down of the fastest rotators in open clusters, 
especially when one reduces the vertical viscosity to a value more or less compatible with the 
surface Li abundance evolution in solar twins. In Section 13. 2^ we present our new dynamic 
model of magnetic braking with the anisotropic turbulent diffusion [v^ ^> u v ) and discuss 
its advantages compared to the original model. Section 13.31 describes a number of stationary 
solutions that address the problem of the penetration of the latitudinal differential rotation 
from the bottom of convective envelope into the radiative interior in the present-day Sun. 
Finally, we summarize our conclusions in Section HI 



2. Simple Models of Rotational Evolution of Solar- Type Stars 

The simple models of rotational evolution of solar- type stars — i.e., the double- zone 
model and th e purely diffusive o ne — have been discussed in detail and compared with 



each other by iDenissenkov et al.l (120091 ). Here, we will only briefly describe the models and 
summarize the main results obtained by comparing their predictions with distributions of an- 
gular velocity of the (assumed) rigidly rotating convective envelope for solar analogs in open 
clusters, Q e = 2n/P, which are derived from the corresponding observed rotation period 
(P) distributions. We will also show that the purely diffusive model fails to simultaneously 
comply with the observational constraints imposed by the spin-down of solar analogs in open 
clusters and the time evolution of the surface Li abundance in solar twins. Although the 
simple models do not identify the physics behind the internal transport of angular momen- 
tum in solar-type stars, they nevertheless provide us with useful estimates of characteristic 
timescales of relevant processes and their dependence on the rotation rate. Therefore, they 
can serve as a good starting point for our analysis. 
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2.1. The Pre-Main Sequence Disk Locking 



In Fig. [U crosses with vertical error bars represent the upper 90th and lower 10th 



percentiles of the fl e distributions for stars having masses in the interval 0.9 < M/M & < 
l.f sampled from open clusters with ages between 2 Myr and 600 Myr (for details, see 
Denissenkov et al.ll2009l ). The data for the youngest clusters set up the initial Q e distributions 
for rotational evolution computations. Of course, these are not the true initial distributions 
inherited by the stars during their formation epoch because they have already been modified 
(reduced) by a disk-locking process compared to what they would have been if the stars' pre- 
MS evolution had taken place in complete isolation. The fact is that after its birth a protostar 
continues to contract before it settles down on the zero-age MS (ZAMS) at an age of about 
40 Myr (for M « 1 M ). During this period of time, the star's radius R and total moment 
of inertia I are decreasing, hence its Q e should be growing as a consequence of angular 
momentum conservation. Contrary to this, it is the surface angular velocity rather than the 
angular momentum that is ob served to be nearly preserved in contracting pre-MS stars (e.g., 
Rebull. Wolff, fc Stroml 12004 ). This is usually explained by a complex magnetic interaction 
that pumps o ut angular momentum from the proto star to its surrounding protoplanetary 
dust disk (e.g. JShu et al.lll994t iMatt fc Pudritzll2005l ) . As a result, the magnetic disk locking 
keeps fl e nearly constant (this is modeled by horizontal fragments of dotted curves in Fig. [Q 
for a period of time up to 10-20 Myr at most before the disk disappears. After that, the 
star's residual contraction and conservation of angular momentum cause an increase of its 
fl P that continues until the star reaches the ZAMS. 



2.2. The Surface Loss of Angular Momentum 

From the moment when the young Sun arrives at the ZAMS (tzAMS ~ 40 Myr) until the 
Sun's present-day age of t@ = 4.57 Gyr, its internal structure does not change significantly. 
In particular, the dimensionless moment of inertia of the Sun's convective envelope i e = 
I e /(R@M@) w 0.0105 remains constant to within 6%. Note that i e amounts only to 14% 
of the Sun's total moment of inertia i = I/(R@Mq) w 0.074. If there were no angular 
momentum transfer from the Sun's radiative core to its convective envelope then the only 
process involved in the Sun's rotational evolution during its life on the MS (izAMS < t < £©) 
would be the braking of envelope rotation by the magnetized solar wind. The corresponding 
rate of angular momentum loss is often approximated by the following equation: 
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where K w is a constant calibrated to give O e (t ) = Q = 2.86 x 10 6 rads 1 (for P & = 
25.4 days), and r2 sat , is a mass-dependent velocity above which the wind gets saturated 



| Chabover. Demarque. fc Pinsonneaultlll995l ; iKrishnamurthi et al.lll997t lAndronov. Pinsonneault. fc Sills 



20031 ). The value of fi sat = 8f2 has been adjusted so that the upper solid curve in Fig. [T] 
simulating the solid-body rotational evolution of rapidly rotating open-cluster solar-type 
stars could fit the uppe r 90th percentiles of their fl e distributions as closely as possible 



(IDenissenkov et al.ll2009l ). 



The differential equation (JTJ can easily be solved. If the initial angular velocity fi 6i zAMS = 
^e(^zAMs) exceeds the saturation threshold fl sat then the solution consists of two different 
parts. The first one, which is valid for tzAMS < t < t aat , is 

^ c = ^c.zams exp [-(t - t ZAMS )/r w ] , (2) 

where the wind braking timescale is r w = (t - £zams)/ (-4/2 + B), while t BaX = t ZAMS + Bt w 
is the time when fl e has decreased to a value of f2 sat . To shorten the equations, we have used 
the notations A = (f2 sat /f2 ) 2 — 1 and B = In (Q e zAMs/ Usat)- The exponential decay of Q c 



is replaced by the Skumanich relation (jSkumanichlll972l ) 



Q = Q& (3) 

y/l + 2 [(A/2 + B)/(l + A)} (t - t )/(t - tZAMS) 

for the time interval t sat < t < t . Finally, the solar-calibrated wind constant has to be 
equal to 

Kw = ^0 M Q (All ± E ) ^ (4) 

Jl| (1 + A) t Q — tzAMS 

Alternatively, if the star arrives at the ZAMS with f2 e ,zAMS < ^sat then its surface angular 
velocity will be declining according to the Skumanich law from the very beginning. The 
above equations can still be used provided that one substitutes f2 sat = fi c ,zAMS m au of 
them. 

Equations ([2]) ^ © suffice to describe the Sun's spin-down for the two limiting cases 
in which its core and envelope rotation are either completely decoupled or, alternatively, 
rigidly coupled. For the first case, which is only of academic interest because it is in conflict 
with uniform rotation in the solar interior, equation (j3J) gives K w = 4.60 x 10 46 cm 2 gs for 
^e,zAMS — 100 Q & , and K w = 4.13 x 10 46 cm 2 gs for fi e ,zAMS = 4.7 Q & (the upper and lower 
solid curves in Fig. [T]). In the second case, which appears to adequately describe the fastest 
rotators, we have to substitute the total moment of inertia i instead of i e into (T4j), resulting in 
K w = 3.24 x 10 47 cm 2 gs and K w = 2.91 x 10 47 cm 2 gs for the same initial angular velocities. 
It is interesting that the solutions ([2]) - ([3]) do not depend on the moment of inertia; therefore, 
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the curves representing them for the two limiting cases completely coincide in Fig. [TJ The 
entire difference between the cases is contained in the values of the wind constant K w . 
These values have to be larger, in proportion to the ratio i/i e = 0.074/0.0105 ~ 7.05, for 
the solid-body rotational evolution because in this case the wind has to be faster to remove 
the correspondingly increased amount of angular momentum, so that finally we could get 
fl e (t&) = ^o- Note that the last constraint comes directly from the observed rotation period 
distributions for solar-type stars in open clusters that show a convergence of P to P Q by the 
solar age. It is also clear that the correct mechanism of angular momentum redistribution 
in solar-type stars should not tolerate a large variation of the wind constant K w between its 
applications to the fastest and slowest rotators as long as we believe that the dependence of 
J on Q e has correctly been factored out in equation (JTJ). 

The illustrated degeneracy of the solutions (J2J) - (EJ) with respect to the two limiting 
cases of coupling between the core and envelope rotation is lifted when we begin to consider 
more realistic cases in which the coupling is allowed to increase progressively in time. To 
demonstrate this via a simple example, let us assume that the moment of inertia of the star's 
rigidly rotating envelope that is decoupled from the rest of it obeys the following rule: 



I(t) 



e.ZAMS 



(5) 



where p = In ( J/J e )/ In (fi^zAMs/^o)- Although this case does not have a physical justifica- 
tion, it admits an analytical solution of equation flTJ and it also guarantees that /(izAMs) — 
(no core/envelope coupling at the beginning) and I(t & ) = I (full coupling at the end). For 
^e.zAMS — 100 f2 , the solution of equation (PJ for the dependence §5§ is plotted in Fig. [T] 
with a dot-dashed curve. It needs K w = 2.56 x 10 47 cm 2 gs. Thus, in spite of the fact that 
the evolution begins with angular momentum being lost only from the convective envelope, 
the wind constant has a value much closer to the previously considered case of solid-body 
rotational evolution. This is explained by the fact that the star will later have to get rid 
of more angular momentum than in the old decoupled case because the moment of inertia 
of its decoupled and uniformly rotating envelope (not just the convective envelope but also 
an outer part of the radiative core adjoint to it) increases well above I e . Therefore, it is 
necessary to lose additional angular momentum from the surface of the convective envelope 
in the future, which will be supplied by the radiative core as its rotation gets more and more 
coupled to that of the envelope, demanding a larger wind constant which, in turn, leads to 
a steep initial decline of Q e . 
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2.3. The Double-Zone Model 



The double-zone model (IMacGregorl Il99ll ) assumes that the Sun (or a solar-type star) 
consists of two uniformly rotating zones, the core and envelope, and that an excess of angular 
momentum in the core, compared to the case in which the whole star rotates rigidly, is 
transferred to the envelope on a specified constant timescale r c . It is evident that, for 
this core-envelope rotational coupling to replenish the angular momentum content of the 
envelope faster than it is drained by the wind, one needs to have r c <C r w = 133 Myr (for 
^e.zAMS = 100), and vice versa. Computations show (e.g.. iDenissenkov et al.ll2009l ) that the 
solid-body rotational evolution of a solar-type star, in the case that its core and envelope 
rotate nearly synchronously all the time, can only be achieved with r c < 1 Myr. 

It turns out that the upper 90th percentiles of Q c distributions (the fastest rotators) 
for solar analogs in open clusters are very well reproduced by the double-zone model with 
r c = 1 Myr (upper dotted curves in panels A and B in Fig. [2]). On the contrary, the slowest 
rotators (the lower 10th percentiles) at ages of order 100 Myr are found to be located below 
the curves representing the solid-body rotational evolution (r c = 1 Myr) that starts from 
the 10th percentiles of Q e distributions for the youngest clusters, even when the safe upper 
limit of Td = 20 Myr is chosen for the disk- locking time (the lower dotted curve in Fig. [2A). 
A rigorous statistical analysis that uses Monte Carlo simulations and compares the full Q e 
distributions (modeled versus observed ones) rather than just their percentiles confirms this 
conclusion. It also gives the following estimates of the coupling time for the slowest rotators: 
r c « 55 ± 25 Myr f or stars with 0.9 <WM Q < 1.1 and t c « 175 ± 25 Myr for stars with 
0.7 < M/M < 0.9 ( iDenissenkov et al.ll2009h . 



2.4. The Purely Diffusive Model 

The purely diffusive model uses appropriate initial and boundary conditions, including 
JT]), to solve the following equation: 

4 dn a ( 4 an\ 

P r ~a7 = ~aZ \ P r v ~aZ I ' ( 6 ) 



dt dr \ dr / 

where v is a constant viscosity whose physical nature is not specified. The time-dependent 
density distribution p(r, t) in the radiative core and other necessary stellar s t ructur e pa- 



rameters are taken from full stellar evolution computations. iDenissenkov et al.l (120091 ) have 
demonstrated that the solid-body rotational evolution of solar-type stars can only be simu- 
lated with v > 10 6 cm 2 s _1 (solid curves in Fig. [2]). They have also established an approximate 
correspondence between the diffusion coefficient v and the coupling time r c from the double- 
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zone model that produce a similar rotational evolution. In particular, this relation says that 
the spin-down of a slowly rotating star computed using the double-zone model with r c w 40 
Myr looks almost identical to that obtained with the purely diffusive angular momentum 
redistribution for z/ ^ 5 x 10 4 cm 2 s _1 (dot-dashed curves in Fig. [2]). This solution gives an 
example similar to (J5]) when the moment of inertia of the rigidly rotating envelope effectively 
increases with time as the diffusive transport of angular momentum gradually penetrates 
into the radiative core. It should be noted that even for a coupling time as long as 90 Myr 
the corresponding diffusion coefficient still remains as large as v « 2.5 x 10 4 cm 2 s -1 . 

The angular momentum redistribution in solar-type stars is known to be accompanied by 
a chemical element transport, the most pronounced manifestation of which is the strongly 



reduced abundance of Li in the solar atmosphere. iMelendez et al.l (120091 ) have recently 
reported that the Sun does not seem to be unique in this respect. It turns out that solar 
twins (stars with mass and metallicity not different from the solar values to within a few 
percent) show a dependence of the surface Li abundance on age to which the Sun also 
belongs (Fig. [3]). A natural interpretation of this relation is that a rate of Li destruction in 
the atmospheres of solar twins is approximately the same, the Sun simply being a relatively 
old star. 

We have added diffusion terms to a network of nuclear kinetics equations relevant for 
the Sun's evolution (the pp chains and CNO cycle reactions). This network has been solved 
in a post-processing way using, as a background, previously stored files with all necessary 
information about the Sun's internal structure from its fully convective pre-MS evolution to 
its present age. The resulting time evolution of the surface Li abundance is plotted in Fig. [3] 
for several values of the diffusion coefficient. It is seen that only a value of v ~ (2-4) x 10 3 
cm 2 s _1 satisfies the observed anti-correlation. On the other hand, a value of v = 5 x 10 4 
cm 2 s" 4 that was used by CHMG93 to model the magnetic braking of solar rotation leads 
to an excessive Li depletion inconsistent with the observations. In our modified model of 
magnetic braking of solar rotation with anisotropic turbulent diffusion, we will employ a 
value of u v = 5 x 10 3 cm 2 s _1 which is close to the values constrained by the Li abundance. 
Note that, without being assisted by another angular momentum transport mechanism, the 
purely diffusive mixing with this small coefficient can reproduce neither the Sun's nearly 
uniform rotation nor the spin-down of the fastest rotators in open clusters (dashed curves in 

Fig. E]). 

Finally, it is important to note that Be, another fragile element, does not seem to share 
the f ate of Li in evo lved MS stars with close to solar effective temperatures, including the 



Sun ( iRandichl 120081 ) . This means that the chemical mixing that causes the Li destruc- 



tion does not penetrate below r w r e — 0.16 R Q ~ 0.55 R Q in the Sun (eg., see Fig. 1 
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of iBarnes. Charbonneau fc MacGregorl Il999l ). r e w 0.713 i? being the radius of its core- 
envelope interface, where the temperature is high enough for protons to begin destroying 
Be. Alternatively, it is possible that th e diffusion coefficient rapidly declines with depth 
(e.g., iRicher. Michaud. &: Turcottdl2000l ). If this is true then the Li mixing may have noth- 



ing to do with the global angular momentum redistribution in the Sun. On the other hand, 
to be in accord with observations the latter process cannot be assisted or accompanied by a 
mass transfer whose rate, when described in terms of a diffusion coefficient, greatly exceeds 



v 



4 x 10 3 cm 2 s 1 . 



3. Dynamic and Stationary Models of Solar Rotation with Large-Scale 

Magnetic Fields 

3.1. Magnetic Braking of Solar Rotation: the Old Solution 



Because our main goal is to modify the well-known model of magnetic braking of solar 
internal differential rotation for it to enable to explain the recent observational data on the 
fl e distributions for solar-type stars in open clusters, without causing a conflict with the 
Li abundance data for solar twins, we will describe the model employed by CHMG93 that 
we want to develop only briefly. For more details about the model, the interested reader 
is referred to the original paper for an excellent analysis of numerical results obtained with 
it and also a discussion of its shortcomings. To save space, we introduce, at the outset, 
the general forms of relevant equations that encompass both the original case considered 
by CHMG93 and our modified case with anisotropic turbulent diffusion. In fact, we have 
chosen an equivalent sy s tem o f equations similar to that used in the follow-up study by 
Ruediger fc Kitchatinovi (119961 ) . which allows the poloidal field potential to fade with time. 



In the spherical polar coordinates {r, 6, if}, the magnetic breaking of solar rotation in 
the presence of anisotropic turbulence in the radiative core is described by the following 
equations 



p r sin 9 



dQ 
~dt 



sin 9 d 

r 3 dr 



pr v v 



dQ 



+ 



P 9 



+ 



1 



Anr 2 sin 9 



dr j rsm 
1 dA d(Br) I 
r 86 



2 9 09 



sin 3 9 



dQ 
~89 



dA d(B sin 9) 



dB 



ld_ 

r dr 



d(Br) 



dr 

Vh_d_ 

r 2 89 



1 

+ - 

r 



dr 

dQ dA dVL dA 
dr d9 d9 dr 



sin 9 dr d9 



1 d(yB sin I 
sin 9 d9 



+ (p r sin 9) SjQ, (7) 
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dA d 2 A sinO d ( 1 dA s 

~dt ~ riv ^ + rih ~^d9\^n~9^9^ (9) 



where Sj = J / (I e £l e ) is a normalized angular momentum loss rate, with I e and Q e being, 
as before, the moment of inertia and angular velocity of the convective envelope. For the 
rate of angular momentum loss fr om the surface J , we co ntinue to use equation ([1]) , whereas 



CHMG93 used the formulation of IWeber & Davisl (119671 ). 



Equations ([7]) — represent the ip - (azimuthal) components of the momentum and 
induction equations. Following CHMG93, we neglect any meridional motion of either ro- 
tational or magnetic origin. We also omit the Coriolis force, while solving the problem 
in an inertial (non-rotating) frame of reference in which the polar axis is directed along 
the rotation axis. The right-hand side of the momentum equation (j7|) contains the viscous 
term, which has been split into its vertical and horizontal components, the magnetic tension 
part of the Lorentz force (l/47r)[(B, V)B] ¥ ,, where the total (poloidal plus toroidal) field is 
B = Bp + Be^ in which the poloidal field has been expressed via its potential A as follows 

„ f Ae v \ f 1 dA 1 dA 1 

B b = Vx — ^ = t-j ¥i ~E Tj-.OL (10) 

\rsmu J {r z sm0 ov rsmv or J 

and, finally, the term simulating the surface angular momentum loss by smearing J over the 
entire convective envelope, which therefore must vanish in the radiative core. 

Equation fl9] ) is decoupled from the fi r st two equations, therefore it can be solved sepa- 



rately. Following iRuediger fc Kitchatinovl (119961 ) . we expand the poloidal field potential into 



a series of the associated Legendre polynomials and keep only the longest living dipole term 

A(r,9,t) » -a{r) sin 2 0exp(-t/r), (11) 

which we substitute into (J9]). As a result, we have to solve the eigenvalue problem 

d 2a n Vh a 
ar r z r 

for the boundary conditions a(0) = a(r b ) = 0, where < r b < Rq, and the potential A will 
thus be determined. 

CHMG93 considered the isotropic case in which v = v y = i\ and rj = i] v = rfc. They 
artificially amplified the microscopic kinematic viscosity v and magnetic diffusivity rj by 
factors of order 10 3 and 10 2 , respectively, making them as large as v = 5 x 10 4 cm 2 s _1 and 
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1] = 5 x 10 5 cm 2 s 1 and keeping their values constant throughout the radiative core. It is 
convenient to transform equation ffl2|) into its dimensionless form 



2 2 ya = -\a, (13) 



d 2 a 2 
dy 2 y 

where y = r/r h , 7 = r/ h /^ v , and A = (rl/rrj v ). Taking the values of 7 = 1 and r/ v = 
5 x 10 5 cm 2 s _1 used by CHMG93, we solve the eigenvalue problem (TT21 [T3|) and find that, 
in this case, the poloidal field potential can be approximated by the following relation: 

2 / \ 2 



A(r,9,t)faB '—[l-—] sin 2 6 x exp(-t/r), (14) 
2 V r bJ 

where r = 35 Myr for 77, = r e = 0.713 -fi^. The same, or very si milar, form of the potential 



was used by both CHMG93 and lRuediger fe Kitchatinovl (Il997l ). However, in both of those 



investigations the time dependence of A was ignored, which might not be a good assumption, 
given the short decay time for A obtained with that for the enhanced magnetic diffusivity. 

The decay time can be as high as 5.4 Gyr, which is comparable to the solar age, only 
when the maximum value of the microscopic magnetic diffusivity in the Sun's radiative core 
77ma.Y = v(t p ) ~ 3 x 10 3 cm 2 s~ 1 is substituted into the above expression for A instead of rj v 



(cf. iRuediger fc Kitchatinovl Il996l ). But, in this case, magnetic braking will hardly work, 
especially, after our having correspondingly decreased the amplified macroscopic viscosity 
either to the maximum microscopic value of z/ max = z/(r e ) w 20 cm 2 s -1 or to its minimum 
turbulent value of v ~ r\ = ?7 max (assuming that r\ is still of turbulent origin but has also 
reached its minimum possible value). The latter substitution brings the turbulent magnetic 
Prandtl number P m = r\jv close to its most plausible value of unity. 

Keeping this inconsistency in mind, we proceed with our revision of the original model 
of CHMG93. Like them, we apply equations ((7j)-(jHl) and ffl~4j) to the whole star assuming 
that v and 77 jump to their large convective turbulent values of order 10 12 cm 2 s _1 at the 
core-envelope interface. Solutions for the angular velocity Q(r,9,t) and toroidal magnetic 
field B(r,9,t) are sought in the space domain (r, 6>) = [0, R Q ] x [0, it] on the time interval 
^zams < t < t Q . Following CHMG93, it is convenient to use /1 = cos 9 as the second 
independent variable instead of 9, in which case the relevant space domain turns into the 
rectangular area (x, fi) = [0, 1] x [—1, 1], where x = r/R Q . 

We start our computations with the following initial conditions: Q(r, 9, 0) = i\zAMS, 
and B{r, 9, 0) = 0. The model's internal structure is considered to be fixed and described by 
the density distribution p{r, t^) for the present-d ay solar model taken from our full stellar 



evolution computations (jDenissenkov et al.ll2009i ). For the dimensionless moment of inertia 



of the convective envelope, we use the time-averaged value of i e = 0.0105. Given that the 
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surface loss of angular momentum is already included in equation ([7]) and that the equations 
are solved for the whole star, the boundary conditions become simple: dQ/dr(R Q ,6,t) = 0, 
dVt/dd = at each point on the rotation axis, and B = on all the boundaries, including the 
core-envelope interface, the latter condition being imposed because of the extremely large 
turbulent magnetic diffusivi ty in the convective envelope t hat makes it feel like a vacuum to 



the internal magnetic field (IRuediger fc Kitchatinovill996l ). 



We consider two poloidal field configurations of the four cases investigated by CHMG93. 
They are specified by the parameter that defines the radial extent of non-vanishing po- 
tential A with respect to r e (equation [TJJ. We have chosen the qualitatively different con- 
figurations D3 and D2 in the notations of CHMG 93. In the most promising, according to 



CHMG93 and iMacGregor fc Charbonneaul (119991 ) . configuration D3 the potential A occu- 



pies the whole radiative core up to the bottom of convective envelope, i.e. we have = r e 
(Fig. [5]A.). The configuration D2 is set up by using the value of = 0.8 R Q > r c , which 
means that A can now penetrate into the lower part of the convective envelope (Fig. [5B) 
but, what is more important, the poloidal field has now a non- vanishing radial component 
at r = r e for some 6. Given that B p in the configuration D2 is thus directly coupled to the 
envelope rotation, it can more effectively transfer into the radiative core the torque applied 
to the convective envelope by the magnetized wind. In the D3 case, there is a time delay 
before envelope braking will be felt in the core. This time is taken for the viscous trans- 
port of angular momentum to build up a sufficiently strong rotational shear just beneath 
the core-envelope interface which will then begin to interact with B p . For a comparison, 
we also introduce a purely diffusive configuration DO with A = (Fig. [7J\). Finally, our 
modified magnetic braking models with anisotropic turbulent diffusion that have the D3 and 
D2 poloidal fields and strong horizontal components of turbulent viscosity in the radiative 
core will be referred to as D3H and D2H. 

To numerically solve the partial differential equations ([7]) - OH]) in the 2D space domain 
with the potential A given by equation f fl^|) we make use of the COMSOL Multiphysics 
software package. Like the custom-made computer code employed by CHMG93, the COM- 
SOL Multiphysics solves PDEs using the finite element method. In our computations, we 
use quadrilateral mesh elements to partition the space domain into a total number of ele- 
ments that slightly varies from one solution to another but that is approximately equal to 
100 x 200 and comparable to the mesh sizes used by CHMG93. Note that we seek a solution 
in a half meridian, whereas CHMG93 did their computations for a quadrant. This reduces 
the resolution of our mesh by a factor of two compared to theirs. Otherwise, apart from the 
different prescription for the surface loss of angular momentum and the option of taking into 
account the exponential decay of B p , our analysis of magnetic braking of solar rotation with 
the isotropic viscosity and magnetic diffusivity is very similar to that by CHMG93. 
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Analyzing the results of their extensive computations, CHMG93 have come to the con- 
clusion that the rotational evolution of the Sun weakly depends on the poloidal field ampli- 
tude B (equations [TD] and [H]). On the contrary, the average strength of the toroidal field 
that is generated by differential rotation in the radiative core changes significantly with B , a 
weaker poloidal field resulting in a stronger toroidal field. This anti-correlation emerges from 
the necessity of the magnetic torque to compensate for the one exerted by the magnetized 
solar wind, the former being proportional to the products -B r -B^ and BqB^ in which B r cx B 
and Bq cx B because the total magnetic field in the radiative core is B = B p + Be^. 

CHMG93 have distinguished the following three important epochs in the magnetic brak- 
ing of solar rotation that differ from each other by the relative strength of the magnetic and 
wind torques. During the first epoch, the magnetic torque is much weaker than the wind 
torque and, as a consequence, the rotational shear produced by the spin-down of the con- 
vective envelope at its interface with the radiative core is spread by the viscosity almost 
without hindrance to the core surface. So, the shear in the outer core grows almost linearly 
with time and, as a result, the toroidal field is amplified there almost quadratically in time. 
This short epoch, lasting between a few thousand to a few million years, is followed by the 
second epoch during which time the differential rotation profile begins to interact with the 
growing toroidal field. For the poloidal field configurations in which the core and envelope 
are magnetically coupled (e.g., the D2 case), this epoch is characterized by vigorous phase 
mixing throughout the radiative core. This process is initiated by large scale toroidal field 
oscillations that originate close to the bottom of convective envelope and propagate through 
the radiative core along the poloidal field lines using them as elastic rubber strings. Because 
the strings are rooted in the convective envelope (for the D2 configuration) at different colat- 
itudes and therefore have different lengths, the Alfven waves carrying the perturbations back 
and forward along the strings get out of phase as time progresses. Eventually, this causes 
large gradients of toroidal magnetic field to build up locally at different places and, as a 
result, the large-scale toroidal field oscillations are diffused away either via oh mic dissipatio n 
or via turbulent magnetic reconnection, depending on the nature of rj (e.g., ISpruit Il999h . 



Finally, during the third epoch, a state of dynamical balance that is achieved at the end of 
the second epoch between the magnetic torque in the radiative core and the total stresses 
(magnetic plus viscous) at the bottom of convective envelope is maintained as the surface 
wind torque decreases with time. 

First, we have repeated the computations of CHMG93 for one of their most favorite 
models, namely the one that starts with f2 ejZ AMS = 50 Q & and has the D3 configuration of 
time-independent poloidal field with the amplitude B = 1 G. Our results are presented in 
Fig. S] (thick solid curves in panels A and C as well as panels B and D). The geometry and 
magnitude of the residual toroidal field in the present-day Sun (panel D) look very similar 
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to those obtained by CHMG93 (upper quadrant of the lower right panel in their Fig. 6). 
As they reported, the model has a nearly flat rotation profile by the solar age (panels B 
and C). However, it seems unlikely that it will be able to fit the upper 90th percentiles of 
the Q e distributions for open-cluster solar analogs after its initial rotation rate is increased 
to a value of fi ej zAMS = 100 Qq corresponding to the observed maximum (the peak on the 
dotted curve representing the solid-body rotational evolution simulated using the double- 
zone model). This supposition is confirmed by our additional test computations that started 
with fi e ,zAMS — 100 o (thin solid curve in panel A). 

It is also not clear what could possibly force the rotational evolution of the CHMG93 
model to change its morphology in the fl e vs. Age plane from the convex one for the fastest 
rotators to the concave one for the slowest rotators, as is indicated by both direct observations 
(e.g., compare the solid and dot-d ashed curves in panels A and B in Fig. [2J respectively) and 



their rigorous statistical analysis (jPenissenkov et al.l 120091 ). Physically, the change implies 



a transition from solid-body rotational evolution to evolution with a differential rotation 
sustained for a few tens to a hundred Myr. As has been mentioned, this cannot be produced 
by a variation of Bq. We can see now that an increase of the initial rotation rate does not 
make the rotational evolution of the CHMG93 model approach the dotted curve either. In 
other words, this model does not have enough internal degrees of freedom to explain the whole 
spectrum of observed rotation rates. Besides, its high viscosity [y = 5 x 10 4 cm 2 s _1 ) causes 
the model to conflict with the observed Li abundances in solar twins (Section 12.41) . When 
we reduce the viscosity to the value of v = 5 x 10 3 cm 2 s -1 , which is more or less compatible 
with the Li data, then the model's rotational evolution becomes absolutely incapable of 
reproducing the upper 90th percentiles of the Q c distributions (dashed curve in Fig. HA). 

As far as the importance of taking into account the exponential decay of the poloidal 
field is concerned, it depends on how big the values of 7 = f]h/f]v and r] v are in the model in 
question. For the above considered case of ^ 6iZ ams — 100 Q Q , 7 = 1 and r] v = 5 x 10 5 cm 2 s -1 , 
a more consistent solution with A decaying on the timescale of r = 35 Myr (equation 
[T41) is plotted in Fig. HA as a dot-dashed curve. Although it does not appear to strongly 
deviate from the old solution (the thin solid curve), the difference between the two solutions 
will probably be significant for Monte Carlo simulations and a rigorous statistical analysis 
of rotation period distrib utions for solar analogs in open clusters like those performed by 



Denissenkov et al.l (120091 ). 



Finally, it is worth noting that the purely diffusive model with the same viscosity v = 
5 x 10 4 cm 2 s -1 as the one used by CHMG93 seems to produce both the Sun's rotational 
evolution and its final internal rotation profile (the dot-dashed curves in panels A and C 
in Fig. [2]) that are not very different from those computed with the model of CHMG93, 
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especially when the poloidal field potential is allowed to decrease with time (the dot-dashed 
curves in panels A and C in Fig. H]). The secondary role of magnetic braking compared to 
the viscous transport of angular momentum is emphasized when we take a much smaller 
value of v = 5 x 10 3 cm 2 s _1 , consistent with the solar Li data (compare dashed curves in 
the same four panels). In this case, the only noticeable effect produced by the presence of 
magnetic fields is a slightly slower rotation of the radiative core at the solar age, especially 
close to its center, which still strongly deviates from the helioseismology data. 



3.2. Magnetic Braking With Anisotropic Turbulent Diffusion 



We start with the aforementioned failed D3 model of magnetic braking of solar rotation 
that has an isotropic viscosity v = 5 x 10 3 cm 2 s _1 and magnetic diffusivity rj = 5 x 10 4 cm 2 s _1 
and which begins its rotational MS evolution with fi e ,zAMS = 100 Q Q (thin dashed curves 
in Fig. E]). Note that its viscosity is already increased by nearly two orders of magnitude 
compared to the maximum value of the microscopic (molecular for the Sun) viscosity in the 
radiative core. This is not a problem though, given that vertical diffusion with a coefficient 
of a similar or slightly stronger magnitude is indeed predicted to occur in ra diative zones of 
solar - metallicity MS stars with M = 1.5 M , especially close to the ZAMS (IPalacios et al. 



20031 ). where internal gradients of the mean molecular weight are too weak to effectively 
hinder it. This diffusion is associated with either rotation-driven meridional circulation or 
small-scale turbulence caused by rotational shear instabilities. It has been successfully used 
to explain chemic al composition anomalies in the atmospheres of massive MS stars and their 
descendants (e.g.. iMeynet fc Maederll2000l ; iDenissenkovl 120051 ). 



Whereas the model of magnetic angular momentum transport based on the hypothetical 
Tayler-Spruit dynamo ha s a problem reproducing differential rotation in the slowest rotators 
(jDenissenkov et al.l 120091 ). the D3 model of CHMG93, on the contrary, cannot provide a 
rapidly rotating solar-type star with a means to spin-down its radiative core fast enough for 
the star's rotational evolution to fit the upper 90th percentiles of the Q e distributions for solar 
analogs in open clusters. The question is whether we can make a reasonable modification of 
the CHMG93 model such that it will become capable of reproducing the rotational evolution 
of both the fastest and slowest rotators while being consistent with the other observational 
constraints? We have shown that a variation of the isotropic viscosity does not help to 
solve the problem self-consistently with or without a magnetic field. The major obstacles 
are the observational limit of v < 5 x 10 3 cm 2 s -1 imposed by the rate of the evolutionary 
Li depletion in solar twins (Fig. [3]) and the rapid decay of the poloidal field caused by high 
viscosity and magnetic diffusivity. On the theoretical side, stellar hydrodynamics does not 
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predict a build-up of vertical turbulent viscosities as large as v ~ 10 6 cm 2 s _1 in the radiative 
cores of rapidly spinning MS solar-type stars, as required for the stars to evolve like solid- 
body rotators. Besides, such large diffusion coefficients would solve the problem without 
being assisted by magnetic fields. 

It should be noted that the poloidal magnetic field does help the viscous transport to 
reduce, even a little, the degree of differential rotation in the Sun's radiative core, as claimed 
by CHMG93. Indeed, panels C in our Figs. [7] and indicate that, whereas the contours of 
constant angular velocity in the present-day Sun's DO model are concentric spheres, they 
become consecutive layers in meridian cross sec tions of an axisymmet ric torus in the D3 



model. This torus, called the "dead zone" by iMestel h Weissl ( I1987I ). is shaped by the 
poloidal field, whose lines form a similar toroidal structure (Fig. EA), according to Ferraro's 
law (Bp, V)f2 = 0. Thus, in the D3 configuration, the viscous transport has to redistribute 
angular momentum throughout a smaller volume than in the absence of magnetic fields (the 
DO case). However, when v = 5 x 10 3 cm 2 s _1 , it does not have enough time to finish its 
work by the solar age. As CHMG93 reported, the situation is slightly improved when the 
poloidal field configuration is switched to D2. In this case, there is no waiting time before 
the spin-down of the convective envelope gets communicated to the radiative core because 
poloidal field lines are now rooted in the envelope, hence they can start transferring the wind 
torque to the core from the very beginning. Besides, because not every poloidal field line 
is now parallel to the bottom of convective envelope (some of them emanate from it, being 
directed inside the core), a larger cylindrical region around the rotation axis (not just the 
axis and its immediate adjacent layers as it was in the D3 case) acquires a rotation rate close 
to that of the envelope. As a result, the dead zone becomes thinner (Fig. [5p). However, 
even in this case the rotational evolution is still far from being close to that of a solid-body 
rotator (thick dashed curves in Fig. [6]). 

When analyzing the above results, especially contour plots C and D in Fig. El we have 
noticed that a strong horizontal turbulence, a concept with different approximate descrip- 
tions that have already been used in other stellar evolution computations, might solve the 
problem of magnetic braking of solar rotation as well. In fact, the same hydrodynamic mod- 
els that give an estimate of v y ~ 10 2 - 10 4 cm 2 s _1 for the vertical component of turbulent 
diffusion also predict a much stronger value, of order 10 4 - 10 6 cm 2 s _1 , for its horizontal 
component in the radiative core of th e solar-metallicity MS star with M = 1.5 M & and 



^e,zAMS ~ 50^0 (IPalacios et all 120031 ) . In a solar- mass star the corresponding turbulent 
viscosities are expected to be lower by a factor of 10, which is roughly proportional to the 
ratio of luminosities for MS stars with M = 1.5 M & and M = M & , because the baroclinic 
and secular Kelvin-Helmholtz instabilities driving this anisotropic turbulence have growth 
rates inversely proportional to the thermal timescale. On the other hand, the coefficient 
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of horizontal turbulent diffusion has been suggested to take on much greater values (by 



more than four orders of magni tude) than previously thought (cf. IZahnlll992l ; iMaederl 12003 



Mathis. Palacios. fc Zahnl 120071 ) 



Zahnl (119921 ) was the first to bring anisotropic turbulent diffusion to the attention of the 



stellar astrophysics community. Anisotropic turbulence has both quite a natural physica l 



justification and an experimental basis in geophysics and oceanography. iTassoull (120001 ) 
notes that "in the Earth's lower atmosphere one has v^jv^ < 10 2 , whereas this ratio may 
be as large as 10 5 in the surface layer of the ocean where large-scale currents are observed". 
It is even more likely to occur in stellar radiative zones with their extremely large Reynolds 
numbers and strong thermal stratifications in the vertical direction. 

Unfortunately, neither v y nor can be derived from first principles yet. Therefore, 
crude estimates are usually adopted and then tuned up by making an order of magnitude 
evaluation of the relevant energy balance, using appropriate results from laboratory or nu- 
merical experiments, and, last but not least, trying to comply with available observational 
constraints concerning mixing in stars. Given this uncertainty, while pursuing our modest 
goal of producing a modified CHMG93 model that should be as simple as the original one, we 
have chosen the following prescriptions for the anisotropic diffusion coefficients in solar-type 
MS stars: 

v v = 5 x 10 3 cm 2 s-\ and u h = 10 6 cm^-i. (15) 

y 100 i ^0 J 

The first expression represents the approximate upper limit for the diffusion rate associated 
with the radial mass transport that is more or less consistent with Li data in solar twins. 
The more realistic value would be closer to 2 x 10 3 cm 2 s _1 (dot-dashed curve in Fig. [3]). 
However, the observational data do not exclude the possibility that u v could be as large as 
3 x 10 4 cm 2 s -1 during the first 300 Myr and then reduced to 1.75 x 10 3 cm 2 s _1 (dotted curve 
in Fig. [3j however, our rotation evolution computations show that this variable viscosity does 
not solve the problem, as is illustrated by Fig. [8]). In other words, our assumed constant 
value for v y should instead be considered a time average, therefore it can be a bit larger, 
say like = 4 x 10 3 cm 2 s _1 (short-dashed curve in Fig. [3]). But it definitely cannot be 
as large as u v = 5 x 10 4 cm 2 s _1 (solid curve in Fig. [3]) and remain constant, as CHMG93 
assumed. As for our choice of i\, its amplitude seems to have a correct order of magnitude 
(in agreement with other most recent applications) for a ZAMS solar model rotating with 



the an gular velocity J\zams — 100 f2 Q , taking into account that iMathis. Palacios. fc Zahn 



(120071) have estimated u h « 10 7 -10 8 cm 2 s 1 for a young (252 Myr) solar- metallicity 1.5 M & 
model star that started its MS evolution with approximately one half the rotation rate. On 
the other hand, the assumed dependence of z/ h on fi e ,zAMS does not seem unnatural and it 
provides us with the necessary additional degree of freedom, such that a variation of its value 
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can probably be used to model the transition from the fastest to slowest rotators, without 
causing a conflict with other observational constraints. 

Thick solid curves in Fig. E] have been computed using the anisotropic diffusion coef- 
ficients (jT5l) for the D2 poloidal field configuration. We have assumed that r] v = 2.5 x 10 3 
cm 2 s _1 and 1% = v\. The poloidal field decay has been taken into account (equation [T4")) with 
an appropriate value of the decay time obtained as a solution of equation ( fT3l) . We can see 
that now the surface rotation evolution curve in panel A fits the upper 90th percentiles of the 
fl e distributions for solar analogs in open clusters, while its corresponding internal rotation 
profile at the solar age appears to be nearly uniform (panel B). A comparison of four similar 
models in panel A demonstrates that it is the combination of the D2 poloidal field geometry 
and strong horizontal turbulence (the D2H model) that results in the rotational evolution 
resembling that of the fastest rotators because neither the D2 nor the D3H models are able 
to produce a convex Q e evolutionary curve. The contours of constant angular velocity and 
toroidal magnetic field for the D2H model of the solar age are shown in panels B and D in 
Fig. [7J It is important to note that, in the absence of internal large-scale magnetic fields, 
strong horizontal turbulent diffusion could not erase the differential rotation in the radiative 
core because it is the dipole poloidal field and its interaction with the rotational shear that 
change the geometry of differential rotation from the spherically symmetric one (Fig. [7_P) 
to the toroidal one (Fig. [5p). On the other hand, the assumption of anisotropic turbulence 
with z/h ^> u y allows us to choose a sufficiently small value of u v , thus avoiding a conflict 
with the observed surface Li abundances in solar twins. 

Finally, when we apply the D2H configuration to model the rotational evolution of solar- 
type stars with initial angular velocities ranging from fl e} zAMS = 100 fi to ^zams = 5 Q , 
we indeed obtain a nice transition from the convex evolutionary curves for the fastest rotators 
to the concave ones for the slowest rotators (Fig. [H]), as hinted by observations. Note that all 
of these models end up with nearly flat internal rotation profiles, and their solar-calibrated 
wind constants do not differ by very much: whereas the fastest rotating model has K w = 
3.64 x 10 47 cm 2 gs, the slowest rotating model needs K w = 4.44 x 10 47 cm 2 gs. 



3.3. Stationary Solutions: the Width of the Solar Tachocline 



The conclusion made by iMacGregor fc Charbonneaul (119991 ) that it is the D3 configu- 
ration, rather than the D2 one, that seems to be the most appropriate one for the model 
of magnetic braking of solar rotation is based on their consideration of the thinness of 
the solar tachocline. Because of the direct magnetic coupling of the convective envelope 
with the radiative core through the D2 poloidal field, the latitudinal differential rotation 
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at the core-envelope interface fl(r e , 9, t ) penetrates into the bulk of the radiative core 
more easily in the D2 case than in the D3 case. The differential rotation of the convec- 



rotating fluid (e.g., 


Kitchatinov & Ruedieer 


1993 


Ruedieer & Kitchatinov 


1996) 


Charbonneau, Dikpati, & Gilman 


( 


1999b and 


Ruediger & Kitchatinov 


( 


2007 


)• 



simple approximate relation 



Q Q (1 -0.15 cos 2 



(16) 



Its values at 9 = 30°, 45°, 60°, and 90° are plotted in panels C and D in Fig. [2] as horizontal 
line segments at r/R > 0.713. It is seen that the variation of the interface angular velocity 
with the colatitude in the present-day Sun is much smaller than the decrease of Q e during the 
Sun's possible previous MS evolution (Fig. [2JA.). Therefore, we have neglected this variation 
in our rotational evolution computations. However, now, when we consider the final model 
of the solar age, we have to be sure that our assumptions and, in particular, our preference 
for the D2 poloidal field configuration, do not lead to the formation of a thick tachocline. 

A stationary solution that has to be tested for the spread of the interfacial differential 
rotation can easily be obtained with the same COMSOL Multiphysics code that we have 
used to compute the rotational evolution of solar-type stars. We have only to set to zero 
the left-hand sides of equations (j7]) - (jHi), omit the term proportional to Sj, use the relation 
ffTB"j) as a new outer boundary condition, and shrink the computational space domain from 



Rq to r 



0.713 R e . We have used this approach to obtain a number of stationary 



solutions, the results of which are discussed below. 

Let us start with non-magnetic solar models. Fig. [TUI A shows a stationary solution 
(contours of constant angular velocity) for one such model in which the kinematic viscosity 
is isotropic and has its near maximum microscopic value of v = 20 cm 2 s _1 that is as- 
sumed (as before) to be constant throughout the radiative core. This solution is evidently 
in conflict with helioseismology data indicating that both the latitudinal and radial differ- 
ential rotation of the solar convective envelope are smoothed out when one crosses the thin 
( Ar < 0.04 R Q ) tac hocline, so that the radiative c ore rotates like a solid body at least down 
to r ~ 0.2 R Q (e.g.. iRuediger fc Kitchatinovll2007l . and references therein). This result is not 
surprising because the isotropic viscosity is known to spread a horizontal velocity inhomo- 
geneity of a given length scale to a radial velocity inhomogeneity of the same length scale 
(IRuediger fc Kit chat inovl 1 19971 ) . Fig. HOB presents a stationary solution for the same non- 
magnetic model but with the anisotropic viscosity whose horizontal component is = 10 4 
cm 2 s _1 ^> u v = 20 cm 2 s" 1 . In this case, the strong horizontal turbulence smooths out the 
interfacial differential rotat ion in a sufficiently t hin layer located immediately beneath the 
core-envelope interface (cf. ISpiegel &: Zahnlll992l ). 



-22 - 



Now, we return to the isotropic viscosity but add a poloidal magnetic field. Fig. ITUD 
shows a stationary solution for the case similar to that considered by CHMG93, namely 
for the D3 poloidal field configuration of amplitude B Q = 0.1 G interacting with the ro- 
tational shear via the (turbulent) viscosity u = 5 x 10 4 cm 2 s _1 , the induced toroidal field 
being diffused away at the rate rj = 5 x 10 5 cm 2 s -1 . In this case, our contours of con- 
stant angula r velocity indeed resemble v e ry clo sely those plotted in the lower right panel 



in Fig. 1 by iMacGregor fc Charbonneaul (119991 ). After returning to the isotropic viscosity 



we would expect the interfacial differential rotation to again be spread radially deep into 
the core but this has not happened. The reason is that now it is the latitudinal trans- 
port of angular momentum by the magnetic stress, proportional to the product BgB v , that 
suppresses the differential rota tion and restricts its penetration into the radiative interior 



(IRuediger fc Kitchatinovl 119971 ). However, the original CHMG93 model did not take into 
account the decay of the poloidal field in spite of the fact that its assumed extremely high 
value of the magnetic diffusivity implies the very short decay time of 35 Myr. With this 
e-folding time, the poloidal field will already have the amplitude of merely 10~ 10 G by the 
age of 725 Myr. Fig. [TOD shows contours of constant angular velocity in the D3 model of the 
solar age with this weak field. Because there is actually no poloidal field left, the solution 
looks almost identical to that for the pure DO case. 

We have chosen the D2 configuration because it more effectively transfers the wind 
torque and rotational shear into the radiative core than the D3 configuration. However, 
this same property plays a negative role when we try to prevent the interfacial differential 
rotation in the present-day Sun from penetrating the radiative interior. Fortunately, our 
D2H model combines the D2 geometry with the assumption of strong horizontal turbulent 
diffusion which, as we have seen in Fig. [TUB, can potentially counteract the negative effect 
produced by the non- vanishing radial component of the D2 poloidal field at the core-envelope 
interface. To test this possibility, we have computed four stationary solutions for our favorite 
D2H model. Fig. [TT1 A shows contours of constant angular velocity at the solar age for 
a D2H model in which neither the anisotropic diffusion coefficients nor the poloidal field 
amplitude were allowed to change with time. This is the worst possible solution in terms 
of the tachocline width. However, like the original CHMG93 solution, our first D2H test 
solution is not fully consistent because it does not take into account the poloidal field decay. 
For our assumed values of r] v = 2.5 x 10 3 cm 2 s _1 and r/h = = 10 6 cm 2 s _1 (for the fastest 
solar-type rotator with f2 e ,zAMS = 100 Q Q ) equation ([TBI gives r = 1.1 Gyr. Hence, by 
the solar age, the initial amplitude B = 0.1 G will be reduced to 2 x 10~ 3 G. Fig. [TTB 
demonstrates that for a D2 poloidal field of such an order of magnitude the tachocline can 
be very thin, due to the continuing strong horizontal turbulent diffusion. Now, given our 
assumption about the dependence of i\ on the angular velocity (equation [15]), it would be fair 
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to construct a solution for the case of having been reduced to 10 4 cm 2 s _1 , proportionally 
to the decrease of fl e from 100 f2 to Q & . Its respective contours of constant Q(r,6,t Q ) are 
presented in Fig.lTTC We see that even in this case the obtained tachocline width is not more 
discrepant than in the original CHMG93 solution (Fig. HUD). Finally, Fig. I11U illustrates 
that the previous solution is greatly improved when the amplitude of the residual poloidal 
field is reduced by another order of magnitude. So, the anisotropic turbulence with the 
horizontal components of turbulent viscosity strongly dominating over those in the vertical 
direction plus the account of the poloidal field decay incorporated into our D2H model make 
it possible to obtain a sufficiently thin tachocline in the present-day Sun's model even for 
the D2 poloidal field configuration. 



4. Conclusion 

The main goal of this work has been to modify the model of magnetic braking of solar 
rotation presented and discussed in detail by CHMG93 so that it could reproduce the most 
recent observational data concerning the rotation of solar analogs in open clusters of different 
ages, without coming into conflict with other available observational constraints. Although 
we have used a different prescription for the surface angular momentum loss rate (equa- 
tion[T]), its corresponding dependence of the envelope braking rate —d(Q e /Q Q )/dt (yr _1 ) on 
the rotation rate fl e /Q Q does not deviate significantly from the dependence predicted by 
the Weber-Davis MHD wind model that was employed by CHMG93 (for comparison, we 



used Fig. 1 from ICharbonneaul Il992l ) . It turns out that even in its original formulation the 
CHMG93 model experiences difficulties in fitting the upper 90th percentiles of the fl e dis- 
tributions for cluster solar-type stars. The discrepancy with the observationally constrained 
rotational evolution of the fastest rotators becomes much worse when we reduce the kine- 
matic viscosity by one order of magnitude and/or take into account the poloidal field decay 
caused by the high magnetic diffusivity in the CHMG93 model. The viscosity reduction 
makes the model more consistent with the data for the atmospheric Li depletion in solar 
twins. The value of v = 5 x 10 4 cm 2 s _1 used by CHMG93 leads to excessive Li destruction 
unless it describes, in terms of an effective diffusion coefficient, another angular momentum 
transport mechanism implemented as a process without a radial mass transfer, e.g. angu- 
lar momentum redistribution by internal gravity waves. However, in this case the question 
would arise as to why that additional process could not do the entire work alone? Poloidal 
field decay also seems to be necessary, given that the model already follows the decay of the 
induced toroidal field. After these reasonable modifications are made, the CHMG93 model 
produces a rotational evolution of a rapidly rotating solar-type star that differs greatly from 
that of the fastest rotators among solar analogs in open clusters. Besides, at the solar age, 
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the model still possesses strong differential rotation in the radiative core, in conflict with 
helioseismology measurements. This differential rotation has such a geometry that its cor- 
responding contours of constant angular velocity form con secutive layers in cross sections of 



an axisymmetric torus, a region called the "dead zone" by lMestel fa Weissl (119871 ) . Its shape 
is a manifestation of Ferraro's law of isorotation, (B p , V)fi = 0, with the poloidal field B p 
having a dipole structure. 

To save the CHMG93 model, we have assumed that the rotation-driven turbulence that 
is thought to amplify the viscosity and magnetic diffusivity compared to their microscopic 
values in stellar radiative zones is actually anisotropic with the horizontal components of 
the transport coefficients strongly dominating over those in the vertical direction. T his is 



not a new hypothesis. It was first introduced to the stellar astrophysics community by IZahn 



(119921 ) and, since then, it has widely been used to model both chemical mixing and angular 
momentum transport in stars. The strong horizontal turbulent diffusion serves a two-fold 
goal in our modified model: it erases the dead zone along isobaric surfaces (a horizontal 
erosion of latitudinal differential rotation), while allowing us to choose a sufficiently small 
value for the coefficient of vertical diffusion that does not lead to a conflict with the Li 
data in solar twins. We have also found it necessary to switch from the D3 poloidal field 
configuration, that was considered to be the most appropriate one in the original CHMG93 
model, to the D2 configuration. In the D2 geometry, the dipole poloidal field partially 
penetrates the convective envelope, which allows it to more effectively transfer the surface 
wind torque to the radiative interior. As a result, a thinner dead zone develops, which makes 
it easier for the horizontal turbulence to erase it. The horizontal component of the turbulent 
viscosity gives us an additional degree of freedom whose assumed linear dependence on the 
initial (ZAMS) angular velocity enables the models to properly reproduce the transition from 
the convex morphology of the surface rotation evolution of the fastest rotators to the concave 
fl e (t) curves for the slowest rotators, as suggested by observations. 



MacGregor fa Charbonneaul (Il999l ) have chosen the D3 geometry as the most suitable 
one for the modeling of magnetic breaking of solar rotation because its corresponding poloidal 
field configuration does not lead to a penetration of the differential rotation of the Sun's 
convective envelope into its radiative interior, consistent with the observations. The helio- 
seismology measurements show that the envelope differential rotation gets smoothed out in a 
thin (Ar < 0.04 R & ) transition layer (the solar tachocline) immediately beneath the bottom 
of the Sun's convective envelope and that the radiative core below the tachocline rotates 
like a solid body at least down to r ~ 0.2 R Q . Contrary to this, it has been shown that 
the non-vanishing radial component of the D2 poloidal field at the core-envelope interface 
causes the interfacial latitudinal differential rotation to move into the bulk of the radiative 
core. However, after our modifications of the original CHMG93 model, the new D2H model 
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appears to be free of this problem. Both the strong horizontal turbulent diffusion and the 
poloidal field decay work towards confining the width of the solar tachocline. 

It is important to note that our D2H model is very simplistic and it is based on a 
number of assumptions whose validity has not been proven. For instance, it is not clear how 
the rotation-driven turbulence should interact with the strong large-scale magnetic fields. 
Also, the formation of the solar tachocline is a much more complex process than described 
by our model. In particular, it must take into account a pene tration of the meridional 



circulation from the conyectiye zone into th e radiative core (e.g., iGough fc Mclntyrdll998 



Ruediger &: Kitchatinovl 120071 ; iGoughl 120091 ) , which our model completely ignores. We are 
aware of these shortcomings. Therefore, we consider our model only as a legitimate extension 
of the CHMG93 model in the sense that we have not made any modifications of the latter 
that are very different from the original assumptions and approximations used by CHMG93. 
On the other hand, we have shown that just one additional assumption, that of strong 
horizontal viscosity, helps to solve several problems with the old model. 
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Fig. 1. — Spin-down of the Sun's convective envelope by the torque from a magnetized 
wind (equation [TJ for the two limiting cases in which the envelope rotation is either totally 
decoupled or, alternatively, rigidly coupled to that of Sun's radiative core (solid curves) 
and for the case in which the moment of inertia of the uniformly rotating envelope, now 
consisting of the convective zone plus a fraction of its adjacent radiative zone, increases with 
time according to the law <^ (dot-dashed curve). Dotted curves (partially overlapped by 
the solid curves) represent the Sun's solid-body rotational evolution computed using the 
double-zone model with the coupling time r c = 1 Myr. Their initial horizontal fragments 
correspond to the pre-MS disk-locking phase. In this and the following plots, the modeled 
evolution of the angular velocity of the Sun's convective envelope is compared with the upper 
90th and lower 10th percentiles of the Q c distributions fo r solar analogs in open cl usters of 
different ages (crosses with vertical errorbars) taken from iDenissenkov et al.l ( 120091 ). 
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Fig. 2. — Spin-down of the surface solar rotation for the case in which the angular momen- 
tum redistribution in the Sun's interior is modeled as a purely diffusive process (equation [5]) 
with a diffusion coefficient v (displayed in panels A and B for the fastest and slowest rota- 
tors, respectively). The corresponding internal rotation profiles at the solar age are shown 
in panels C and D with crosses r eprese nting the observational data from helioseismology 
measurements by ICouvidat et al.l (j2003l ). Horizontal line segments at r/R > 0.713 in the 
same panels show the latitudinal differential rotation of the (bottom of) Sun's convective 
envelope at 9 = 30°, 45°, 60°, and 90° (from the lower to upper). 
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Fig. 3. — Evolutionary decline of the atmospheric Li ab undance in solar twins. Symbols 
represent the observational data of iMelendez et al.l ( 120091 ). Curves are our solar model post- 
processing computations of the kinetics of nucleosynthesis in the pp chains and CNO cycle 
with chemical mixing modeled as a purely diffusive process. For Li to start being destroyed, 
it has to be mixed at least down to r « 0.65 R Q (in the present-day Sun, the core-envelope 
interface is located at r = r e 0.713 Rq). For mixing penetrating below r w 0.6 R Q , the 
surface Li depletion depends almost entirely on the value of the diffusion coefficient u and 
not on the mixing depth (curves). The dotted curve represents the possibility of v being 
enhanced during the first 300 Myr. 
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Fig. 4. — Magnetic braking of solar rotation computed using the D3 poloidal field config- 
uration (Fig. [5]A.) with the isotropic viscosity v = 5 x 10 4 cm 2 s _1 and magnetic diffusivity 
7] = 5 x 10 5 cm 2 s" 1 (thick solid curves in panels A and C as well as contours of constant 
angular velocity and toroidal field in panels B and D, the last three panels showing solutions 
at the solar age). These values are equal to those used in the favorite model of CHMG93. 
In the models represented by dot-dashed, thin solid and dashed curves, the initial angular 
velocity was increased from 50 Q & to 100 fl , and in the first of them the poloidal field of the 
same amplitude B Q = 1 G as in CHMG93 was allowed to decay with time (equations [T3"1- [T4"|) . 
The CHMG93-type solution obtained with the reduced values of rj = 5 x 10 4 cm 2 s _1 and 
v = 5 x 10 3 cm 2 s _1 , the latter being more or less compatible with the observed Li depletion 
in solar twins, fails to reproduce both the spin-down of the fastest rotating solar-type stars 
in open clusters and the internal rotation profile of the present-day Sun (dashed curves in 
panels A and C). 
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Fig. 5. — Geometries of the poloidal field potential, that is assumed to have here the constant 
amplitude B = 0.1 G, in the configurations D3 (panel A) and D2 (panel B), and their 
corresponding final (solar age) contours of constant angular velocity and toroidal field (panels 
C and E, and D and F, respectively) obtained in our computations of magnetic braking of 
solar rotation with the isotropic values of v = 5 x 10 3 cm 2 s _1 and rj = 5 x 10 4 cm 2 s _1 . Note 
the strong residual differential rotation in the torus-shaped "dead zones" (panels C and D). 
Red circles in panels A and B show the core-envelope interface. 
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Fig. 6. — Comparison of four similar cases of magnetic braking of solar rotation for the 
poloidal field configurations D3 and D2 with (suffix "H") and without a strong horizontal 
turbulent diffusion (u^ ^> z/ v ). In the cases D2H and D3H, z/ v = 5 x 10 3 cm 2 s _1 , r] v = 2.5 x 
10 3 cm 2 s -1 , and the poloidal field of the amplitude B = 0.1 G is allowed to decay. In the D2 
and D3 cases, we used v = 5 x 10 3 cm 2 s _1 , rj = 5 x 10 4 cm 2 s _1 , and a constant poloidal field 
of the same amplitude. In the D2H and D3H cases, we assumed that u h = r] h = 10 6 cm 2 s _1 . 
Panel B compares the final internal rotation profiles with the helioseismology data. Note 
that only our D2H model gives a satisfactory solution for the fastest rotators (thick solid 
curves) . 
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Fig. 7. — Panel A: a sketch of the DO model of a solar-type star without magnetic fields, 
the red circle representing the bottom of its convective envelope. Panels B and C: contours 
of constant angular velocity in our D2H model and in the non-magnetic model at the solar 
age. Panel D: contours of constant toroidal magnetic field in the final D2H model. 
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Fig. 8. — Illustration that the assumption of v being enhanced during the first 300 Myr 
(dotted curve in Fig. [3J does not help to solve the problem using the D3 poloidal field 
configuration. 
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Fig. 9. — Solid curves show the rotational evolution of solar-type stars computed with 
our D2H model using the anisotropic turbulent viscosity ( Fl5l) . Other model parameters 
are t] v = 2.5 x 10 3 cm 2 s _1 (except for the slowest rotating model that has r] v = 5 x 10 3 
cm 2 s _1 ), ?7h = z/h, B = 0.1 G initially and is decaying exponentially with an e-folding time 
that is calculated using equation ffl3|) . The solar-calibrated wind constant varies between 
K w = 3.64 x 10 47 cm 2 gs and K w = 4.44 x 10 47 cm 2 gs from the fastest to slowest rotator. 
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Fig. 10. — Penetration of the latitudinal differentia l rotation of the bottom of S un's con- 
vective envelope, Q(r e , 9, t )/fi = (1 — 0.15 cos 2 6) (IRuediger fc Kitchatinovl 120071 ). into its 
radiative core. DO: v = 20 cm 2 s -1 has its near maximum microscopic value. D0H: same 
as DO but with the horizontal components of = 10 4 cm 2 s" 1 , note the formation of a 
thin tachocline in this case, as originally proposed by ISpiegel fc Zahnl (Il992l ). D3CHMG93, 
Bq = 0.1 (G): the original D3 model of CHMG93 with the c onstant poloidal field (cf. lower 
right panel in Fig. 1 from iMacGregor fc GharbonneaiJfl999h . D3CHMG93, B = lO" 10 (G): 
same as in the previous panel but with a strongly decayed poloidal field amplitude, which 
had not been taken into account by CHMG93. 



-39 - 




Fig. 11. — Same as in Fig. dO] but for our magnetic braking model (the D2 poloidal field 
with the anisotropic turbulent diffusion). D2HSTAT, B = 0.1 (G): this case assumes that 
all the model parameters remain constant, including the value of = 10 6 cm 2 s _1 (for 
^e,zAMS = 100 Q Q in equation IT5]) . D2HSTAT, B = 0.001 (G): same as in the previous 
panel but with the decayed poloidal field. D2H, B = 0.001 (G) and B = 0.0001 (G): here, 
it is additionally assumed that u h has been reduced to 10 4 cm 2 s _1 , proportionally to the 
decrease of fl e from 100 Q G to Q Q . 



